Identifying the Magnetoelectric Modes of Multiferroic BiFeOs 
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We have identified three of the four magnetoelectric modes of multiferroic BiFeOs measured using 
THz spectroscopy. Excellent agreement with the observed peaks is obtained by including the effects 
of easy-axis anisotropy along the direction of the electric polarization. By distorting the cycloidal 
spin state, anisotropy splits the ^/±i mode into peaks at 20 and 21.5 cm -1 and activates the lower 
<3>±2 mode at 27 cm -1 (T = 200 K). An electromagnon is identified with the upper mode at 
21.5 cm -1 . Our results also explain recent Raman and inelastic neutron-scattering measurements. 

PACS numbers: 75.25.-j, 75.30.Ds, 78.30.-j, 75.50.Ee 
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Multiferroic materials hold tremendous technological 
promise due to the coupling between the electric po- 
larization and cycloidal magnetic order. Two classes of 
multiferroic materials have been identified. In "proper" 
multiferroics, the cycloid develops at a lower tempera- 
ture than the ferroelectric polarization; in "improper" 
multiferroics, the electric polarization is directly coupled 
to the cycloid []], Q so they develop at the same tem- 
perature. Although the multiferroic coupling is typically 
stronger in "improper" multiferroics, those materials also 
have rather low transition temperatures. Due to the high 
magnetic transition temperatures of "proper" multifer- 
roics like BiFe0 3 with T N « 640 K [3-5], those materials 
continue to attract a great deal of interest. With the 
availability of single crystals, the physical understanding 
of BiFeOs has advanced rapidly over the past few years 
but the description of the distorted spin cycloid and mi- 
croscopic interactions in BiFeOs remains incomplete. 

Recent inelastic neutron-scattering measurements on 
single-crystal samples of BiFeOs S 0] were used to es- 
timate the easy-axis anisotropy K along the polariza- 
tion direction as well as the Dzyaloshinskii-Moriya (DM) 
coupling D that induces the long-period cycloidal or- 
der with wavevector Q = (27r/a)(0.5 + (5,0.5,0.5 — 5) 
and S ~ 0.0045 (8l4lo|. However, more precise estimates 
for those coupling constants may be obtained based on 
the magnetoelectric peaks recently measured using THz 
spectroscopy [11]. In this work, we evaluate the mag- 
netoelectric peaks from a model that includes easy-axis 
anisotropy and DM interactions. In addition to provid- 
ing excellent agreement for three of the four experimental 
peaks, we also provide several new predictions for the se- 
lection rules governing those three peaks and one higher- 
energy peak. 

Due to the displacement of the Bi 3+ ions, ferroelectric- 
ity appears in rhombohedral BiFeOs below T c « 1100 K 
12j . For the pseudo-cubic unit cell in Fig.l with lat- 




FIG. 1: (Color online) The pseudo-cubic unit cell for BiFeOs 
showing x', y ; , and z directions. The distorted spin spiral 
propagates along the x' direction with spins in the (— 1, 2, — 1) 
plane. AF interactions Ji and J2 are also indicated. 



tice constant a w 3. 96 A, cycloidal order develops be- 
low Tn with propagation vector (27r/a)(5, 0, —S) along x' 
for each (1, 1, 1) plane [3|-[5j. Neighboring (1, 1, 1) planes 
are coupled by the antiferromagnetic (AF) interaction 
Ji between the S = 5/2 Fe 3+ spins. The DM interac- 
tion D along y' or [—1,2,-1] produces a cycloid with 
spins in the ( — 1,2,-1) plane. Easy-axis anisotropy K 
along [1,1,1], parallel to the polarization P, distorts the 
cycloid by producing odd harmonics of the fundamental 
ordering wavevector [l3j |. 

Although Tn is much lower than T c , the magnetic do- 
main distribution of BiFeOs can be effectively manipu- 
lated by an electric field 0, 0, EH . In a magnetic field 
above 20 T, the transformation of the cycloid to an al- 
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most commensurate structure with weak ferromagnetic 
moment [15j is accompanied by a sharp drop of about 40 
nC/cm 2 in the electric polarization 13, 16, I?}- There- 
fore, the additional polarization P md observed below Tn 
is induced by the cycloidal spin state. 

This coupling between the cycloid and electric polar- 
ization is produced by the inverse DM mechanism [181-120] 
with induced polarization P md oc e^- x (Si x Sj) where 
eij = Rj — Ri and S^ are the Fe 3+ spins. Within each 
(1,1,1) plane, e^- = \^2ax.' connects spins at R^ and 
R^ = Hi + V^ax'. So Si x Sj points along y' and P ind 
points along z'. 

To evaluate the spin- wave (SW) excitations of BiFeOs, 
we use the Hamiltonian 0, Qj], [2l[ : 

H = j 1 j2s i -s j + j 2 j2 s i- s j - ^E( s * • z ') 2 
-^E w+v , ax ,y'-^x^)- (D 

In the first and second exchange terms, denotes a 

sum over nearest neighbors and (i, j)' a sum over next- 
nearest neighbors. The third term originates from the 
easy-axis anisotropy along z' and the fourth term from 
the DM interaction with D along y'. For a fixed cy- 
cloidal period, D is a smoothly increasing function of K, 
as shown in Fig. 5 of Ref.[7[. 

Since 6 = 0.0045 is close to 1/222, a unit cell con- 
taining 222 sites within each (1,1,1) plane was used to 
characterize the distorted cycloid, which was expanded 
in odd harmonics of the fundamental wavevector Q = 
(27r/a)(0.5 + £,0.5, 0.5 -£) 
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S Z ,(R) = S^2 CWi cos ((2m + 1)Q • R), (2) 

ra=0 



K = meV 



K = 0.002 meV 



S X >(R) = VS 2 -MR) 2 sgn(sin(Q • R)). (3) 

The harmonic coefficients obey the sum X^m=o ^2m+i = 
1. 

Excitation frequencies and intensities were evaluated 
by performing a 1/S expansion in the rotated frame of 
reference for each spin [23| in the 444-site unit cell with 
two AF-coupled layers. Due to zone folding, there are 222 
positive eigenfrequencies for each wavevector. However, 
only a handful of those frequencies have any intensity. 

To evaluate the SW frequencies, we use AF interactions 
J\ =4.5 meV and J2 = 0.2 meV, which describe the 
inelastic neutron-scattering measurements 0, [?J at 200 
K Q. Below 4 meV and for < q < 25, all possible SW 
frequencies are denoted by the dashed curves in Fig. 2 for 
K = and 0.002 meV. Branches with nonzero intensity 
are indicated by the dark solid curves. 

When K = 0, the two strongest SW branches are 
plotted in the inset to Fig. 2 (a) with points at multi- 
ples of the wavevector q = S labeled as <l> n or \I/ n . For 



> 

CD 




FIG. 2: (Color online) The SW modes of BiFe0 3 versus q/8 
for wavevector (27r/a)(0.5 + q, 0.5, 0.5 — q). Dashed lines show 
all possible excitations and the solid lines show only those 
modes with non-zero intensity, (a) For K = 0, the two 
strongest SW branches are plotted versus q/8 in the inset 
with points <£ n and ^ n at multiples of 8 for c = 1. The solid 
(green) point at ^1 indicates a MR and EM mode, (b) For 
K = 0.002 meV, the \£±i mode (green and blue) splits and 
the lower $±2 mode (red) is activated. The x (blue), y (red), 
and z (green) components of the MR modes are indicated. 
The EM mode corresponds to the upper (blue) \£±i mode 



small frequencies and q = n8, cj(<£ n _i) 
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uj(^ n ) = \A + Ti 2 cJ, where c is the SW velocity of the 
linear branches [25|. In the reduced-zone scheme, odd-n 
\I/ n and even-n <£ n modes lie at the zone center q = S 
while even-n \I/ n and odd-n <I> n modes lie at the zone 
edge q = 0. 

When K > 0, the SW spectrum for small frequencies 
changes dramatically. Higher harmonics of the cycloid 
split every set of crossing &± n and ^± n modes. The 
largest splitting occurs at q = between the $±1 modes. 
A smaller splitting occurs at q = 5 between the \£±i 
modes. While the $0 mode is shifted slightly above zero 
frequency, the ^0 mode is moved up to just below the top 
$±1 mode. Figure 3 plots the evolution of those points 
with anisotropy. Although too small to appear in this 
plot, even the $±2 modes are split by anisotropy. 

Spectroscopic intensities are given by the matrix ele- 
ments of M a = gjiB Sia for magnetic resonance (MR) 
and by the matrix elements of the induced polarization 



)ind 



OC 



{x'x(S,xS,)} a 



(4) 



for the electromagnon (EM) [26|. Each matrix element 
is evaluated between the ground state |0) and an excited 
state \q) containing a single magnon with wavevector at 
the cycloidal zone center q = S. 

When K = 0, there is a single MR peak at the ^1 point 
in Fig.2(a). For this mode, |(g|M^|0)| = |(g|M^|0)| > 
and |(g|M y /|0)| = 0. An EM peak with component y' 
coincides with the MR peak. Only the y' component of 
(g|P^ nd |0) is nonzero. No other frequencies at q = 5 are 
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al. [TlJ conjectured that the 20 and 21.5 cm -1 lines were 
produced by the splitting of the \I/±i modes due to a 
modulated DM interaction along z'. We conclude that 
the splitting of the \I/±i modes is caused by easy-axis 
anisotropy along z' with quantitatively accurate values. 
Perhaps due to the small matrix element plotted in the 
inset to Fig- 3, the lower $±2 peak at 27 cm -1 was not 
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FIG. 3: (Color online) The evolution of the SW frequencies 
with easy-axis anisotropy. Solid and dashed pairs of curves 
indicate split ±n modes. Inset are the matrix elements for the 
MR modes (arbitrary units) versus K. The horizontal dashed 
lines indicate the experimental results of Ref.[ll|] at 200 K. 



magnetoelectrically active for K = 0. 

When K > 0, both \I/±i modes are active with nonzero 
x' and z' MR matrix elements for the upper and lower 
modes, respectively. The MR matrix elements |(g|M a |0)| 
of those modes are plotted versus anisotropy in the inset 
to Fig. 3. While the matrix element of the lower ^±1 
mode decreases with anisotropy, that of the upper ^±1 
mode increases. Because the third harmonic of the \£±i 
modes couples to the first harmonic of the $±2 modes, 
the lower $±2 mode is activated by anisotropy [25] with 
a nonzero y' MR matrix element plotted in the inset to 
Fig.3. 

A large EM peak with matrix element (g|P^ d |0) co- 
incides with the upper mode. Another EM peak 
with matrix element about 10 times smaller was found 
at the upper ^±3 mode with frequency 5.39 meV (43.4 
cm -1 ). This peak also has a significant MR matrix ele- 
ment |(g|Ma;/|0)|. Hence, the active ^±(2m+i) peaks cor- 
respond to out-of-plane modes of the cycloid (excited by 
magnetic fields along x' or z') while the active $±2™ 
peaks correspond to in-plane modes of the cycloid (ex- 
cited by a magnetic field along y'). 

Below 30 cm -1 , THz spectroscopy [11] observed four 
infrared modes with wavevectors (measured at or extrap- 
olated to 200 K 0) of 17.5, 20, 21.5, and 27 cm -1 . The 
upper three mode frequencies are denoted by the dashed 
lines in Fig.3. These modes can be quite accurately de- 
scribed by K = 0.002 meV, which produces MR peaks 
at 2.49, 2.67, and 3.38 meV, remarkably close to the ob- 
served peaks at 2.48, 2.67, and 3.35 meV. Talbayev et 



detected [ll| above about 150 K. 

The selection rules governing the MR modes also agree 
with the results of Ref.[ll]. Field directions hi || [1, — 1, 0] 
and h.2\\ [1, 1, 0] can be written as hi = x'/2 — y/3y' /2 and 
ti2 = x'/2 + V3y' /6 + y // 2/3z / . Consequently, the lower 
$±2 mode with MR component y' and the upper 
mode with MR component x' are excited both by fields 
hi and I12, but the lower \I/±i mode with MR component 
z' is only excited by field h2 [27]. The predicted upper 
\I/±3 mode with MR component x' should be excited by 
both fields hi and Ii2. 

Since it is associated with a large EM peak, the up- 
per mode can also be excited by an electric field 
along y' . The observation of non-reciprocal directional 
dichroism (NDD) under an external magnetic field B ex 
along z' [26] would confirm this prediction. NDD re- 
quires linearly-polarized electromagnetic waves propa- 
gating along x' with electric and magnetic components 
E w ||y' and H w ||z'. 

By contrast, the observed low-energy mode at wavevec- 
tor 17.5 cm -1 (2.17 meV) cannot be explained by our 
model because there are no zone-center excitations be- 
low the lower ^±1 mode at 2.5 meV. Nevertheless, the 
2.17 meV peak lies tantalizingly close to both the up- 
per <£±i mode and the ^0 mode at the zone bound- 
ary in Fig. 2(b). It is possible that the 2.17 meV peak 
is associated with a more complex magnetic structure 
induced by some interaction or anisotropy not included 
in our model. Indeed, additional long-range order with 
wavevector (27r/a)(5, 0, — S) or (27r/a)(0.5, 0.5, 0.5) would 
hybridize the $±1 and ^0 modes at the zone edge with 
the ^±1 modes at the zone center. In Raman spec- 
troscopy, the anomalies of the 2.17 meV peak at 140 K 
and 200 K were identified with spin reorientations 28l.l29| 
of the cycloid. However, any spin reorientation must re- 
spect the selection rules for the zone-center MR modes. 

All of the active zone-center modes have been observed 
with Raman spectroscopy [28l-[30| . The Raman peak [loj 
at 23 cm -1 (T = 10 K) detected with parallel polar- 
izations (exciting out-of-plane modes) can be identified 
with the ^±1 modes and the 25.5 cm -1 peak detected 
with crossed polarizations (exciting in-plane modes) can 
be identified with the lower $±2 mode. Another Raman 
peak [29] at 43.4 cm -1 (T = 80 K) lies quite near the 
estimated wavevector of the upper ^±3 mode. 

Because it is governed by different selection rules than 
THz spectroscopy, Raman spectroscopy has also observed 
several zone-edge cycloidal modes. That includes strong 
peaks at 34.5 cm" 1 (4.28 meV, T = 80 K) [29] and 
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FIG. 4: (Color online) The averaged neutron-scattering inten- 
sity x a t (27r/a)(0.5, 0.5, 0.5) versus uj for values of K from 
to 0.005 meV. Measurements at T = 200 K are in the inset. 

32 cm- 1 (3.96 meV, T = 10 K) [30], close to the pre- 
dicted energy of the zone-edge \I/±2 modes. As conjec- 
tured above, additional long-range order may produce 
such zone-edge cycloidal peaks in the Raman spectrum. 

Earlier estimates for K and D were based on low- 
energy inelastic neutron-scattering measurements taken 
at the wavevector (27r/a)(0.5, 0.5, 0.5) [7]. Because the 
inelastic-scattering cross section is quite broad, estimates 
based on those measurements are not as precise as the 
estimates described above based on spectroscopy mea- 
surements. By using the predicted frequencies of the $±i 
modes to match the primary inelastic peaks at 1.1 and 
2.5 meV, Matsuda et al. estimated that K w 0.005 meV 
[24j . However, peaks in the inelastic intensity x'^cj) may 
be slightly shifted away from those frequencies due to the 
finite resolution of the measurements. 

Averaged over a realistic resolution function and in- 
cluding the three sets of twins states with wavevectors 
along x', y', and z', x"^) ls plotted in Fig. 4 for six val- 
ues of K from to 0.005 meV. The inset provides the 
inelastic measurements at 200 K, which are almost tem- 
perature independent between 100 and 300 K. For K > 0, 
the lowest two peaks are near the 3>±i points. Higher- 
energy peaks might be expected near the ^±1 and $±2 
zone-center modes and near the ^±2 zone-edge modes. 
But due to the instrumental resolution, those peaks are 
smeared into broad features in x" i 00 )- 

The best qualitative agreement with the inelastic mea- 
surements is obtained using K w 0.004 meV. Below 5 
meV, the measured x"^) contains peaks at 1.2, 2.4, 3.4, 
and 4.4 meV. While four nearby peaks are predicted by 
our model, the lowest-energy peak is too weak and the 
two highest-energy peaks are slightly too low compared 
to the experimental results. 

A value for K between 0.002 and 0.004 meV is consis- 



tent with the prediction K = 0.0027 meV obtained from 
Monte-Carlo simulations of the phase diagram in a mag- 
netic field [Uj. As discussed above, anisotropy produces 
higher harmonics C2 m +i>i of the cycloid. Within the 
predicted range of K, 624 > (C1/C3) 2 > 176. Elastic 
neutron-scattering [H and NMR measurements (3l[ indi- 
cate that /1//3 = (C1/C3) 2 is given by 500 and 25, re- 
spectively. Because more precise estimates are provided 
by THz spectroscopy, we believe that a value of K closer 
to the lower limit of 0.002 meV and a value for (C1/C3) 2 
close to 500 is likely. 

We conclude that the upper three peaks measured by 
THz spectroscopy for BiFe03 can be quantitatively pre- 
dicted by a model that includes both DM and easy-axis 
anisotropy interactions. Our results for the MR compo- 
nents in Fig. 3 can be used to check the selection rules 
governing those three peaks. A higher-energy EM peak 
at 5.39 meV is predicted with MR component x' . The 
low-energy 2.17 meV peak measured by THz and Raman 
spectroscopies may be produced by additional long-range 
order. 
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